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Abstract 

The role of the surface during polarization switching in constrained ferroelectrics 
is investigated using the time-dependent Ginzburg-Landau theory. The model 
incorporates the elastic and electrostrictive effects in the form of a long-range 
interaction that is obtained by eliminating the strain fields subject to the elastic 
compatibility constraint. A square shaped finite sized constrained ferroelectric 
system with vanishing surface polarization is considered. Computer simulations of 
the hysteresis reveal that the corners of the constrained square act as nucleation 
sources for 90° domain structures. It is also found that no nucleation of 90° domain 
structure takes place in a system with semi-infinite (corner free) geometry. For the 
corner free case, the polarization switches by 180° reorientations and a front of 
the reverse domains emerges from the surface, eventually sweeping over the entire 
system. Size dependence of the hysteresis loops is also studied. 
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I. INTRODUCTION 



The spontaneous polarization in ferroelectrics can be switched by applying a strong electric 
field opposite to the polarization direction It is crucial to understand the details of this 
switching process due to many applications in memory devices and electromechanical systems. It 
is now a well established fact that the switching behavior is strongly influenced by nucleation, 
growth and motion of domain walls. Due to the nucleation events, the experimentally measured 
coercive field is much smaller than the intrinsic value predicted by the simple 1-D Landau theory 
0. In a recent work ||, we qualitatively showed that dipolar defects can act as nucleation centres 
for 90° domain structures. However, for a constrained system, such as a grain in a ferroelectric 
ceramic, the surface itself is a heterogeniety that can cause localized nucleation. Although, there 
have been some studies investigating the role of a surface on ferroelectric phase transitions [§,§J, 
the polarization switching in finite sized systems has not been well studied theoretically. 

In this paper, we study the influence of a surface during polarization switching using the 
time dependent Ginzburg-Landau model described in |J. In particular, we are interested in 
understanding the role played by the surface in nucleating 90° ferroelastic domains during the 
switching process. Recent experiments on BaTiOs single crystals [f|] have demonstrated the 
importance of 90° domains during polarization reversal. The theory incorporates elastic and 
electrostrictive effects in the form of an anisotropic long-range interaction obtained by eliminating 
strain fields, subject to the elastic compatibility requirement. Instead of using periodic boundary 
conditions (corresponding to a bulk system), we consider a finite sized constrained ferroelectric 
with the boundary condition of vanishing polarization at the surface. Physically, this corresponds 
to a system where the free charge at the surface is compensated. It is believed that this charge 
compensation occurs at the grain boundaries in ferroelectric ceramics. Due to the above discussed 
boundary conditions, we can neglect the effects of the depolarization field. The nonlocal interaction 
is represented by the appropriate polarization gradients terms in the free energy expansion. 

The organization of this paper is as follows. In section 2, we give details of the time dependent 
Ginzburg-Landau model, starting from a free-energy functional for the polarization field. Section 
3 gives details of our simulations of polarization switching for a finite square shaped system. In 
section 4, we describe simulations of switching behavior in a semi-infinite constrained system. We 
end the paper by a summary and discussion of our results in section 5. 

II. MODEL 

The 3-D Ginzburg-Landau free energy for ferroelectric systems has been given in Ref. |7|]. In 
this work, we restrict to a 2-D system undergoing a square to rectangle transition, for which the 
total free-energy is given by 

F = J df\f t + f 9 + f d + f es + f ext ] (1) 
where f] is the local free-energy density given by 

ft = + p y 2 ) + X (P * 4 + PyA) + ^f p * p ^ ( 2 ) 

where P x and P y are the components of the polarization vector, a\, an and a\2 are the free energy 
parameters. The term f g represents the gradient energy, 
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where gi, g2 and g 3 are the gradient parameters. It has been shown that the gradient energy may 
be obtained as a local approximation to the nonlocal interactions |TJ]. The term f e i represents 
the usual elastic energy of the square system. We consider the bulk strain 4>i = (r) xx + r] yy )/\^2, 
deviatoric strain 2 = (Vxx — Vyy)/^ an d shear strain <p 3 = r\ xy = rj yx . Here rjij is the linear 
elastic strain tensor given as rjij = ^(ff 1 + §^r), where Ui, (i = x,y) represents the dispacement 
components. The elastic free-energy can then be written as 



lei = y^i + y02 + —<P3 , 



(4) 



where ai, a 2 and a 3 are constants that can be expressed in terms of linear combinations of elastic 
constants of the system. The coupling energy f es may be written as, 
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[P 2 + P y 2 ) - q 2 <\> 2 {P 2 - Py 2 ) ~ qsfoPxPy 



(5) 



Here qi, q 2 and q 3 are constants that are related to the electrostrictive constants of the system. 
The term f ext is given by f ext = —E ext ■ P. The amplitude of field E ext is a tunable parameter in our 
simulations for studying polarization switching. We assume that the system reaches mechanical 
equilibrium very fast so that we may integrate out the elastic fields, subject to the elastic compat- 
ibility constraint f|. In terms of the strain components 0i, <p 2 and 3 , the elastic compatibility 
relation is given as 
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(6) 



To incorporate the elastic compatibility constraint, we consider an effective elastic part of the free 
energy F e ff = F e i + F es + Peons, with Lagrangian multiplier A, where 



F el = J dk{^\Mk)\ 2 + f \Mk)\ 2 + ^\Mk)\ 2 }, 



(7) 



Fes = - J dk( [ q 1 r 1 (k)M-k) + g 2 r 2 (£)0 2 (-£) + g 3 r 3 (£)0 3 (-£) j (8) 

Fcons = J dk^i-eM-k) + {k 2 - k y 2 )<j) 2 (-k) + V8k x k y <p 3 (-k))}- (9) 

Here k 2 = k 2 + k 2 and Ti(k) , T 2 (k) and T 3 (k) are, respectively, the Fourier transforms of 
P 2 + P y 2 , P 2 — P y 2 , and P x P y . The condition of mechanical equilibrium and the compatibility 
relation demands = 0(i = 1, 2, 3) and = 0. 

oi^i(jfe) = giri(jfe) + k 2 X{k) 
a 2 (j) 2 (k) = q 2 T 2 (k) - (k x 2 - k y 2 )\(k) 

a 3 fo(k) = q 3 T 3 {k) - V8k x k y \{k), (10) 

subject to the constraint 
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-k 2 M-k) + {k 2 - k 2 )<p 2 (-k) + V8k x k y( f) 3 (-k) = 0. (11) 

Using the above equations, the multiplier X(k) is given as 

- _ -fc 2 giri(-fc)/ai + (k 2 - k y 2 )q 2 T 2 (-k)/a 2 + \f&k x k y q 3 T 3 (-k)/a 3 
1 j " k*/a 1 + (k x 2 -ky 2 y/a 2 + 8k x 2 ky 2 /a 3 ' 1 ' 

At this stage we make a simplification by assuming that the polarization electrostrictively couples 
only to the deviatoric strain. This assumption will not effect the qualitative behavior of the model 
in context surface effects, which are the main focus of this paper. We should also keep in mind 
that for the present case of a square to rectangle ferroelectric transition, the deviatoric strains 
are the dominant mode of deformation. The bulk and shear strains are usually localized at the 
domain walls for this case. Including the coupling with bulk and shear deformations may introduce 
intermediate metastable states with more complex patterns. Exclusion of these couplings makes 
the computation faster due to faster equilibriation. Finally, we should also mention that the 
neglect of this coupling in the free energy does not meen that bulk and shear strains will not be 
generated in this model at all. In fact, in this model, the strains are related to each other through 
the elastic compatibility relations which leads to the creation of bulk and shear strains, even under 
the above described approximation. In the limit of no coupling of polarization with shear and 
bulk strains, i.e. q\ — > and q 3 — > 0, using the equilibrium conditions and the lagrange multiplier 
\(k), the free energy F e ff can be expressed in Fourier space as 

2 



F eff = ^-J dkH(k)\T 2 (k)\ 2 , (13) 

where H(k) = + h 2 2 {k) - 2h 2 {k) + ^M). The quantities h^k) = k 2 Q{k), h 2 (k) = {1 - 

{k 2 — k y 2 )Q(k)} and h 3 (k) = —y/8k x k y Q(k), with Q(k) defined as 

m = (kya + (k! 2 -ky 2 yl8k x 2 ky 2 /py (14) 

We have introduced dimensionless constants a = ai/a 2 and f3 = a 3 /a 2 . 

The effective interaction derived above is strongly direction dependent and is crucial to describe 
the domain wall orientations. Similar anisotropic interactions have been considered in context of 
martensitic transformations PJIU]- 



As we wish to study heterogeneous nucleation, a dynamical formalism is needed to describe the 
non-equilibrium effects associated with domain switching. The time dependent Ginzburg-Landau 
model for the polarization fields is used. We introduce rescaled time variables t* = (t\cei\L)(ai < 



0)( L is the kinetic coefficient ), space variable r* = r/9, the wave vector k* = k6 {6 = J g\/a\oc\\), 
where a is a dimensionless constant. The polarization field is transformed as P x = Pru and 



P y = Prv, where Pr = \J\ai\/a±i is the remnant polarization of the homogeneous state without 
elastic effects. With this set of parameters, the dimensionless time dependent Ginzburg-Landau 
equations are given as 

u, t * = u — u 3 — duv 2 + au, x * x * + bu,y* y * + cv, x * y * + e x — 7« / dk*H(k*)T(k*)exp(—ik* ■ r*) 

v, t * = v — v 3 — dvu 2 + av, y * y , + bv, x , x , + cu, x , y * + e y + jv J dk*H(k*)T(k*)exp(—ik* ■ r*). (15) 

Here, e = E ext /(Pn\ai\) is the rescaled external electric field. The notation u, t represents the 
time derivative. The defined constants are: b = (g 2 /\ai\9 2 ), c = (g 3 /\ai\9 2 ), d = (a^/an) and 
7 = (<?2 2 / '"11^2) • The quantity T(k*) is the Fourier transform of {u 2 — v 2 ). 
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III. SIMULATION OF POLARIZATION SWITCHING FOR A CONSTRAINED 

SQUARE SHAPED FERROELECTRIC 



We now describe the details of our simulations on the switching behavior for a finite constrained 
ferroelectric. The long-range interaction is conveniently defined in Fourier space, which allows us 
to use the Fourier Transform method. This enforces the need to have periodic boundary conditions 
at the surfaces. To define a finite system using this scheme, we consider an 88 x 88 square region 
within a full simulation grid of size 128 x 128. The constraining is implemented by assuming 
that the polarization variables u(r*) and v(r*) vanish at the boundary and outside of the 88 x 88 
region. This configuration represents an electromechanically clamped ferroelectric in which all 
displacements vanish at the surface. Equation (15) is discretized using the Euler scheme on 
the above described grid. The space discretization step Ax* = Ay* = 1 and the time interval 
At* = 0.02. The free energy parameters chosed by us are d = 2 and 7 = 0.05. The gradient 
coefficients are chosen asa — b — c — 1. Similarly, we choose the elastic parameter ratios as 
a — f3 — 1. We first start from a random small amplitude initial conditions for the fields u(r*) 
and v(r*) and an external field e x * = 10. We solve equation (15) to obtain a single domain 
configuration with the polarization oriented towards the +x* gradually decaying at the surface of 
the square system. Figure 1(a) shows this configuration corresponding to an electric field value 
e x * = 10. In this figure, the arrows represent the local polarization vector ( we have not shown all 
the polarization vectors for the sake of legibility). Notice the gradual decay of the amplitude of 
the vectors in the surface region. This configuration is used as the initial state for the simulations 
of polarization switching. Equations (15) are now solved at each time step, starting from the 
initial configuration depicted in figure 1(a) with a time dependent external field e x *(t*). The 
time dependence is given as e x * = eo[l — 2t*/T], where eo = 10 is the initial electric field value 
and T = 5000 is the total number of simulated switching time steps. As the field decreases, the 
magnitude of the polarization vectors decreases uniformly every where in the simulated system at 
the beginning (keep in mind that the polarization is fixed to zero at the surface for all times). As 
the external electric field approaches zero, we observe that near two of the corners, the polarization 
vectors are distorted towards the [11] directions. Figure 1(b) shows the pattern corresponding to 
e x * = 0.22 where we can see the diagonal orientation of the polarization vectors near the corners. 
This indicates that shear strains are generated near the corners. The result is in agreement with 
recent simulations by Jacobs for a constrained ferroelastic system undergoing a square to rectangle 
transition where it was demonstrated that nondeviatoric strains are generated at the corners 



11||L2|1 . These diagonal shear strains at the corners are responsible for nucleation of 90° domain 
walls, as seen figure l(c)(e x * = 0.199) and figure l(d)(e x * = 0.0). In figure 1(e) [e x * = —0.199), we 
can clearly see a well developed twinned structure with 90° domain walls. As we further decrease 
the field to e x * = —0.38, domains polarized along [01] direction grow at the expense of domains 
polarized along [10] directions, as can be seen by comparing figure 1(e) and figure 1(f). In figure 
1(f) {&x* = —0.38), we can also observe that the dipoles in the vicinity of the corners and domain 
walls have started reorienting towards the — x* direction. Thereafter, the reversed domains grow 
at the expense of domains polarized along [01], and eventually all the dipoles are oriented along 
the —x* direction. This growth can be clearly seen from figure 1(g) (e x * = —0.46) and figure 1(h) 
(e x * = —10). Figure 1 illustrates the special role played by the corners during switching. The 
corners act as nucleation centres for 90° domain structures. Unlike the classically held belief of 
180° reorientation of polarization, the switching takes place by sucessive 90° reorientations of the 
dipoles. Recently, evidence for such non 180° switching has been observed in rhombohedral single 
crystals of PZN-PT [|13j. The reverse cycle can be simulated by using the configuration in figure 
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1(h) as the initial condition and increasing the field as e x * = eo[l — 2t*/T], where eo = — 10 and 
T = 5000 is the total number of switching time steps. Figure 2 shows the hysteresis loops obtained 
by plotting the average polarization < u > with the field e x *. The solid line corresponds to the 
situation depicted in figure 1 and the dashed line corresponds to this loop for the case without 
any surface (periodic boundary conditions). We can observe that for the case of constrained 
square, the simulated hysteresis loop shows a long "tail" region where the forward and reverse 
cycles curves do not close up. This corresponds to the appearence and disappearence of the 90° 
domain pattern. In contrast, for the periodic boundary condition case, hysteresis loop has no "tail" 
region as there are no nucleation events. We can also observe that the saturation polarization 
and remnant polarization field are lower for the constrained system than that for the case with 
periodic boundary conditions. However, there is not much difference between the coercive fields 
of the two cases. We have simulated the hysteresis loops for different system sizes to study the 
size dependence of the switching phenomena. In figure 3, we show the simulated hysteresis loops 
for system sizes indicated at the top of each frame. A comparison of the loops for sizes L = 88 
and L = 68 reveals that the coercive fields and the saturation polarization are almost the same 
for both systems. Also, the 90° domain structure is shorter lived for the smaller system since the 
"tail" region of the hysteresis loop for the L = 68 system is smaller than that for the L = 88 
system. For smaller sizes, the domain pattern becomes increasingly shortlived during the switching 
process, as can be concluded from the loop for L = 48 where the tail has moved further inwards. 
However, the saturation polarization and the coercive fields still do not show much variation. The 
tail almost disappears for the L = 20 hysteresis loop and we see a small bump near the coercive 
field region. For this size, we can see a reduction in the coercive field as well as the saturation 
polarization due to the surface effects. Below this size, nucleation of 90° domains does not take 
place. The switching occurs by 180° reorientations that preferentially starts at the surface. The 
hysteresis loops in this regime have shape similar to the periodic boundary condition case. It 
is clear from figure 3 that on further decreasing the size, the saturation polarization, remnant 
polarization and coercive field decrease rapidly. We can see that the hysteresis is very small for 
L = 6 which is the smallest size simulated. For the present set of parameter values, we do not 
observe a ferroelectric to paraelectric transition for L = 6. Such a transition was observed for a 
simulation with much higher value of the gradient coefficients. The behavior observed in figure 3 
can be understood in terms of size dependence of ferroelectric behavior. It is experimentally and 
theoretically established that the spontaneous polarization progressively decreases on decreasing 
size for a system in which the polarization is constrained to decay at the surface []5T,|6T,|i~4^|15 . Thus, 



the remnant polarization, saturation polarization and coercive field decrease on decreasing the 
system size. The anisotropic nonlocal interaction responsible for the nucleation of 90° domain 
structures also decreases as the size is decreased. Thus for small sizes, the hysteresis loops are free 
of 90° nucleation. A similar transition from multidomain to single domain states on decreasing the 
size has been discussed in earlier works |l^Jl6[| . The main effect of size reduction comes from the 
surface layer under the vanishing polarization boundary condition. The transition temperature is 
effectively reduced and the energy barrier between the low temperatures states is also lowered. 



IV. SIMULATION OF POLARIZATION SWITCHING FOR A SEMI INFINITE 

CONSTRAINED FERROELECTRIC 

In this section, we describe switching behavior for a semi-infinite constrained ferroelectric 
system. Here we use periodic boundary conditions along the y direction and vanishing polarization 
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conditions at the surface for x direction. Figure 4 shows the dynamics of polarization switching for 
the semi-infinite case. Except the boundary conditions, all other parameter values, including the 
time dependence of the switching field are the same as that for the case of the constrained square 
described in section 3. Figure 4(a) shows the fully polarized single domain initial configuration 
for e x * = 10. In figure 4(b), we can see that a single domain persists even when the applied field 
is lowered to e x = —0.319, although the polarization magnitude has decreased by a very small 
amount. As the field is further decreased to e x * ~ —0.36, nucleation of the reversed polarization 
domains from the two surfaces takes place, as depicted in figure 4(c). The reversed domains 
grow larger and eventually into a single domain crystal with reversed polarization. This reversed 
domain growth is clear from figure 4(d) (e x * = —0.399) and figure 4(e) (e x * = —0.44). Figure 
4(f) shows the final state corresponding to e x * = —10. Figure 4 illustrates the fact that for the 
corner free semi-infinite geometry, the reversal occurs only by 180° reorientations. The hysteresis 
loop corresponding to the situation shown in figure 4 can be computed. In figure 5, we plot the 
hysteresis loop for the case shown in figure 4 along with the loop for the surface free case with the 
periodic boundary conditions. Unlike the case shown in figure 2, the hysteresis loops for the surface 
free case and the semi-infinite case have the same shape. However, the saturation polarization, 
the remnant polarization and the coercive field are slightly smaller for the semi-infinite case due to 
the suppressed transition at the surface. We have also simulated the size dependence for this case 
and found that the hysteresis loops become progressively narrower with decreasing size without 
any shape change. This is due to the fact that there are no 90° domains and the long-range elastic 
interaction does not play any role in this case. 

V. SUMMARY AND DISCUSSION 

We have simulated the switching behavior in constrained ferroelectrics with the boundary 
condition of vanishing polarization at the surface. A square shaped system is considered along 
with a system with semi-infinite geometry. We find that the corners of the square shaped system 
nucleate 90° domains and the switching is dominated by the motion of 90° domain walls. In 
contrast, for the semi-infinite geometry, there is no nucleation of 90° domains and the switching 
proceeds by 180° reversal. Clearly, the corners in the constrained square play a crucial role during 
the switching process. We believe that the anisotropy of the long-range elastic interaction interacts 
with the corners to produce inhomogenieties that result in 90° nucleation. Our results suggest 
that the switching behavior and the hysteresis loops in real systems may depend on the shape of 
the grain and its boundary conditions. 

We have also simulated the size dependence of hysteresis loops. The decrease of the satura- 
tion polarization, remnant polarization and the coercive field with decreasing size is common to 
both configurations. This can be understood in terms of the supression of the transition in an 
electromechanically constrained system as the system size approaches the correlation length of 
order parameter. The free energy for a constrained system is a function of the system size. For 
the constrained square case, the polarization reversal process is dominated by intermediate 90° 
domain pattern changes for larger systems and becomes totally via 180° reorientations for size 
less than a critical dimension. This transition is due to the weakening of the long-range elastic 
interactions with decreasing system size. 

Finally, one should note that the behavior of decreasing coercive field with decreasing size 
indicated by the simulation results does not apply to thin films || , but only to clamped grains in 
ferroelectric polycrystals. For thin films, not all the dimensions are reduced simultaneously, which 
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is different than our condition here. In addition, the misfit strain with the substrate is important, 



which has to be incorporated in the theory [17]. 
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Figure captions: 

Figure. 1: Pattern evolution for the simulated switching in a constrained square. The snap- 
shots correspond to (a) e x * = 10; (b) e x * = 0.22; (c) e x * = 0.199; (d) e x * = 0; (e) e x * = —0.199; 
(f) e x , = -0.38; (g) e x . = -0.46; (h) e x , = -10. 

Figure. 2: Plots of < u > vs e x showing the simulated hysteresis loops for the constrained 
square (solid line) and the periodic boundary condition case (dotted line). 

Figure. 3: Simulated hysteresis loops for different sizes of the constrained square. 

Figure. 4: Pattern evolution for the simulated switching in a semi-infinite constrained system. 
The snapshots correspond to (a) e x * = 10; (b) e x * = —0.319; (c) e x * = —0.36; (d) e x * = —0.399; 
(e) e x . = -0.44; (f) e x . = -10. 

Figure. 5: Plots of < u > vs e x showing the simulated hysteresis loops for the semi-infinite 
constrained system (solid line) and the periodic boundary condition case (dotted line). 
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